from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *
import mydataMN as mydatmn
import Functions as func

def iext1(t):
    if t <= 500*ms:
        return 0*uA*cm**-2
    elif t <= 5500*ms:
        return (t/ms-500) * 25/5000*uA*cm**-2
    elif t <= 100000*ms:
        return (25-(t/ms-5500) * 25/5000)*uA*cm**-2
    else:
        return 0*uA*cm**-2
def i_inj(t):
    md = mydatmn.MNVars()
    return t/second * 0*uA*cm**-2
#    if t <= 100*ms or t>= 1500 *ms:
#        return 0*nA
#    elif t > 1000*ms:
#        return -20*uA*cm**-2 
#    else:
#        return 20*uA*cm**-2 

def current_inj(t):
    if t > 0*msecond:
        return -0.5*namp
    else:
        return 0*mamp
    
        
class NeuronFactory:
    
    def generate(self, size, init1='def', init2='def'):
        md = mydatmn.MNVars()
        md.__init__()
        
        eqsMN=Equations('''
        dvm/dt = (iext1(t)-(md.gna*mna*mna*mna*hna*(vm-md.Ena))-(md.gk*mk*mk*mk*mk*(vm-md.Ek))-ICa-(md.gkca*ca/(ca+md.kd)*(vm-md.Ek))-(md.gl*(vm-md.Elmn))-(md.gc/md.p*(vm-vmd)))/md.cpf : mvolt
        dmna/dt = (func.m_inf_na(vm/mvolt)-mna)/md.taumna : 1
        dhna/dt = (func.h_inf_na(vm/mvolt)-hna)/func.tau_h_na(vm/mvolt) : 1
        dmk/dt = (func.m_inf_k(vm/mvolt)-mk)/func.tau_m_k(vm/mvolt) : 1
        dmcan/dt = (func.m_inf_can(vm/mvolt)-mcan)/md.tau_m_can : 1
        dhcan/dt = (func.h_inf_can(vm/mvolt)-hcan)/md.tau_h_can : 1
        ICa = md.gcan*mcan*mcan*hcan*(vm-md.eCa) : mA*umetre**-2
        dca/dt = md.f*(-md.alpha*(ICa)-md.kca*ca) : mole*dmetre**-3
        
        
        dvmd/dt = (-(md.gnapDe*mnapd*hnapd*(vmd-md.Ena))-ICad-(md.gkcaDe*cad/(cad+md.kd)*(vmd-md.Ek))-(md.glDe*(vmd-md.Elmn))-(md.gcDe/(1-md.p)*(vmd-vm)))/md.cpf :mvolt
        dmnapd/dt = (func.m_inf_nap(vmd/mvolt)-mnapd)/md.tau_m_nap : 1
        dhnapd/dt = (func.h_inf_nap(vmd/mvolt)-hnapd)/func.tau_h_nap(vmd/mvolt) : 1
        dmcand/dt = (func.m_inf_can(vmd/mvolt)-mcand)/md.tau_m_can : 1
        dhcand/dt = (func.h_inf_can(vmd/mvolt)-hcand)/md.tau_h_can : 1
        dmcald/dt = (func.m_inf_cal(vmd/mvolt)-mcald)/md.tau_m_cal : 1
        ICad = (md.gcanDe*mcand*mcand*hcand*(vmd-md.eCa))+(md.gcalDe*mcald*(vmd-md.eCa)) : mA*umetre**-2
        dcad/dt = md.f*(-md.alphaDe*(ICad)-md.kca*cad) : mole*dmetre**-3
        
        ''')
        #eqsMN+=OrnsteinUhlenbeck('I',mu=0.0*nA,sigma=0.025*nA,tau=10*ms)
        MN=NeuronGroup(1,model=eqsMN,
        threshold=EmpiricalThreshold(threshold=-20*mV),
        implicit=True,method='exponential_Euler')
        
        
        # Initialization
        vin =md.vinit1
        if init1 != 'def':
            vin=init1
        MN.vm=vin*mV
        MN.mna = func.m_inf_na(vin)
        MN.hna = func.h_inf_na(vin)
        MN.mk = func.m_inf_k(vin)
        MN.mcan = func.m_inf_can(vin)
        MN.hcan = func.h_inf_can(vin)
        MN.ca = -md.alpha*(md.gcan*MN.mcan*MN.mcan*MN.hcan*(MN.vm-md.eCa))/md.kca
        
        vin = md.vinit2
        if init2 != 'def':
            vin=init2
        MN.vmd=vin*mV
        MN.mnapd = func.m_inf_nap(vin)
        MN.hnapd = func.h_inf_nap(vin)
        MN.mcand = func.m_inf_can(vin)
        MN.hcand = func.h_inf_can(vin)
        MN.mcald = func.m_inf_cal(vin)
        MN.cad = -md.alpha*((md.gcan*MN.mcand*MN.mcand*MN.hcand*(MN.vmd-md.eCa))+(md.gcal*MN.mcald*(MN.vmd-md.eCa)))/md.kca

        return MN